{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing a sparse solution of a set of linear inequalities\n",
    "\n",
    "A derivative work by Judson Wilson, 5/11/2014.<br>\n",
    "Adapted from the CVX example of the same name, by Almir Mutapcic, 2/28/2006.\n",
    "\n",
    "Topic References:\n",
    "\n",
    "* Section 6.2, Boyd & Vandenberghe \"Convex Optimization\" <br>\n",
    "* \"Just relax: Convex programming methods for subset selection and sparse approximation\" by J. A. Tropp\n",
    "\n",
    "\n",
    "## Introduction\n",
    "\n",
    "\n",
    "We consider a set of linear inequalities \n",
    "$Ax \\preceq b$ \n",
    "which are feasible. We apply two heuristics to find a sparse point $x$ that satisfies these inequalities.\n",
    "\n",
    "The (standard) $\\ell_1$-norm heuristic for finding a sparse solution is:\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{minimize}   &  \\|x\\|_1 \\\\\n",
    "    \\mbox{subject to} & Ax \\preceq b.\n",
    "    \\end{array}\n",
    "\n",
    "The log-based heuristic is an iterative method for finding\n",
    "a sparse solution, by finding a local optimal point for the problem:\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{minimize}   &  \\sum_i \\log \\left( \\delta + \\left|x_i\\right| \\right) \\\\\n",
    "    \\mbox{subject to} & Ax \\preceq b,\n",
    "    \\end{array}\n",
    "where $\\delta$ is a small threshold value (which determines if a value is close to zero).\n",
    "We cannot solve this problem since it is a minimization of a concave\n",
    "function and thus it is not a convex problem. However, we can apply\n",
    "a heuristic in which we linearize the objective, solve, and re-iterate.\n",
    "This becomes a weighted $\\ell_1$-norm heuristic:\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{minimize}   &  \\sum_i W_i \\left|x_i\\right| \\\\\n",
    "    \\mbox{subject to} & Ax \\preceq b,\n",
    "    \\end{array}\n",
    "which in each iteration re-adjusts the weights $W_i$ based on the rule:\n",
    "    $$W_i = 1/(\\delta + \\left|x_i\\right|),$$\n",
    "where $\\delta$ is a small threshold value.\n",
    "\n",
    "This algorithm is described in papers:\n",
    "\n",
    "* \"An affine scaling methodology for best basis selection\"<br>\n",
    "  by B. D. Rao and K. Kreutz-Delgado\n",
    "* \"Portfolio optimization with linear and fixed transaction costs\"<br>\n",
    "  by M. S. Lobo, M. Fazel, and S. Boyd\n",
    "\n",
    "## Generate problem data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "\n",
    "# Fix random number generator so we can repeat the experiment.\n",
    "np.random.seed(1)\n",
    "\n",
    "# The threshold value below which we consider an element to be zero.\n",
    "delta = 1e-8\n",
    "\n",
    "# Problem dimensions (m inequalities in n-dimensional space).\n",
    "m = 100\n",
    "n = 50\n",
    "\n",
    "# Construct a feasible set of inequalities.\n",
    "# (This system is feasible for the x0 point.)\n",
    "A  = np.random.randn(m, n)\n",
    "x0 = np.random.randn(n)\n",
    "b  = A.dot(x0) + np.random.random(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\ell_1$-norm heuristic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "status: optimal\n",
      "Found a feasible x in R^50 that has 40 nonzeros.\n",
      "optimal objective value: 28.582394099513873\n"
     ]
    }
   ],
   "source": [
    "# Create variable.\n",
    "x_l1 = cp.Variable(shape=n)\n",
    "\n",
    "# Create constraint.\n",
    "constraints = [A*x_l1 <= b]\n",
    "\n",
    "# Form objective.\n",
    "obj = cp.Minimize(cp.norm(x_l1, 1))\n",
    "\n",
    "# Form and solve problem.\n",
    "prob = cp.Problem(obj, constraints)\n",
    "prob.solve()\n",
    "print(\"status: {}\".format(prob.status))\n",
    "\n",
    "# Number of nonzero elements in the solution (its cardinality or diversity).\n",
    "nnz_l1 = (np.absolute(x_l1.value) > delta).sum()\n",
    "print('Found a feasible x in R^{} that has {} nonzeros.'.format(n, nnz_l1))\n",
    "print(\"optimal objective value: {}\".format(obj.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterative log heuristic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 1: Found a feasible x in R^50 with 48 nonzeros...\n",
      "Iteration 2: Found a feasible x in R^50 with 36 nonzeros...\n",
      "Iteration 3: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 4: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 5: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 6: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 7: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 8: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 9: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 10: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 11: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 12: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 13: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 14: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 15: Found a feasible x in R^50 with 33 nonzeros...\n"
     ]
    }
   ],
   "source": [
    "# Do 15 iterations, allocate variable to hold number of non-zeros\n",
    "# (cardinality of x) for each run.\n",
    "NUM_RUNS = 15\n",
    "nnzs_log = np.array(())\n",
    "\n",
    "# Store W as a positive parameter for simple modification of the problem.\n",
    "W = cp.Parameter(shape=n, nonneg=True); \n",
    "x_log = cp.Variable(shape=n)\n",
    "\n",
    "# Initial weights.\n",
    "W.value = np.ones(n);\n",
    "\n",
    "# Setup the problem.\n",
    "obj = cp.Minimize( W.T*cp.abs(x_log) ) # sum of elementwise product\n",
    "constraints = [A*x_log <= b]\n",
    "prob = cp.Problem(obj, constraints)\n",
    "\n",
    "# Do the iterations of the problem, solving and updating W.\n",
    "for k in range(1, NUM_RUNS+1):\n",
    "    # Solve problem.\n",
    "    # The ECOS solver has known numerical issues with this problem\n",
    "    # so force a different solver.\n",
    "    prob.solve(solver=cp.CVXOPT)\n",
    "    \n",
    "    # Check for error.\n",
    "    if prob.status != cp.OPTIMAL:\n",
    "        raise Exception(\"Solver did not converge!\")\n",
    "\n",
    "    # Display new number of nonzeros in the solution vector.\n",
    "    nnz = (np.absolute(x_log.value) > delta).sum()\n",
    "    nnzs_log = np.append(nnzs_log, nnz);\n",
    "    print('Iteration {}: Found a feasible x in R^{}'\n",
    "          ' with {} nonzeros...'.format(k, n, nnz))\n",
    "\n",
    "    # Adjust the weights elementwise and re-iterate\n",
    "    W.value = np.ones(n)/(delta*np.ones(n) + np.absolute(x_log.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Result plots\n",
    "\n",
    "The following code plots the result of the $\\ell_1$-norm heuristic, as well as the result for each iteration of the log heuristic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAGoCAYAAAATsnHAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3X90XPV55/HPI1kgDLZlCcsJdYg9giAn3WBkiWSh2RZLpu1Jf4KMs02bPe1imbQ93Z5tIkG3Tbv7xzpSON1uTpvGcna77Wl2CxY0aU9PQyUbkrbJJvphkg0BQzwmjZMmNpbH/AoY28/+ce/I49HIGt258tw7836do6OZO3dmnkE+8+He7/c+X3N3AQCQNA3VLgAAgFIIKABAIhFQAIBEIqAAAIlEQAEAEomAAgAkUlUCysxOmdm0mQ0XbOs3sz4zG6xGTQCAZFlRpffd4e4T+Ttm1i9J7j5hZhkz6yt8HABQf6p1iq/FzDIF93skZcPbWUldl78kAECSVOsIqlXSrJntdffdklqKHm8rfoKZDUgakKSrr756a2dn5/JXCQCI3fT09Avuvm6x/aoSUO4+KklmlgtP7+UUhNZizxmVpO7ubp+amlr2OgEA8TOzb5Wz32U/xWdmA/kxJ0knw9+TunAUlZE0frnrAgAkSzXGoB6WlDOzPkly9zF3H5OUKdjGBAkAqHOX/RSfu+ck5QNoomD7yOWuBQCQXFyoCwBIJAIKAJBIBBQAIJEIKABAIhFQAIBEIqAAAIlEQAEAEomAAgAkUioD6o1zXu0SAADLLJUBdeTEy3InpACglqUyoN44d15P/8tL1S4DALCMUhlQknTwme9XuwQAwDJKZUBd1dSoA88cr3YZAIBllMqAWt3cpCe/ndMLL79e7VIAAMsklQG16qoVcpeeOHyi2qUAAJZJKgPqqqZGrV99JeNQAFDDUhlQkrSts11fePYFnTl7vtqlAACWQYoDar1efv2sJp+frXYpAIBlkNqAuv2GNl2xokEHnmY2HwDUotQG1MorVui2jjYdeOb7dJUAgBqU2oCSpN7Odn3r5KvKvvBKtUsBAMQs1QF1R2e7JOkgp/kAoOakOqA2rF2pzjet0gGmmwNAzUl1QEnBdPPJ50/p9A/eqHYpAIAYpT6geje369x51xeepasEANSS1AfUlres1dqVTTpI81gAqCmpD6jGBtMdN7Xr8cPHde48080BoFakPqAkadvmduVefUOH/vlUtUsBAMSkJgLqPTeu04oGY40oAKghNRFQa65qUs/GVq6HAoAaUhMBJQWz+Q5//yV9e/bVapcCAIhBzQTUtrCrxOOHOYoCgFpQMwGVWXeNNl17Nd3NAaBG1ExAScFR1JeOnNQrr5+tdikAgArVVED1drbrzLnz+qdvvlDtUgAAFaqpgOre2KpVV65gHAoAakBNBdQVKxr0b962TgeePs4ihgCQcjUVUFIwDnX8pdf11HdfrHYpAIAK1FxA/dhN62QmZvMBQMrVXEC1XXOlbnlLiw6yiCEApFrNBZQk9W5er68eO63jL71W7VIAABHVZEDlu0o88QyLGAJAWtVkQHW+aZWuW9OsA5zmA4DUqsmAMjNt29yuf3juBb1+9ly1ywEARFCTASVJvZ3r9eqZc/pydrbapQAAIqjZgPrXHW1qbmrQQRYxBIBUqtmAam5q1I/ccK0OPPN9ukoAQArVbEBJ0rbO9fr27A/0zeMvV7sUAMAS1XhABdPND3CaDwBSp6YD6k1rmvWO61brIG2PACB1ajqgpGCNqKlvzSr36plqlwIAWIKaD6htm9frvEuff5auEgCQJjUfUO/8oTW69por6G4OAClT8wHV0GC646Z2PXH4uM6eO1/tcgAAZar5gJKk3s3tevG1s5r+1qlqlwIAKFNdBNSP3LhOTY1GVwkASJG6CKhrrlyhd2fauB4KAFKkLgJKCi7a/ebxl/Wtk69UuxQAQBnqKqAkcZoPAFKibgLqrW1X64b2awgoAEiJugkoKegq8X+zJ/Xy62erXQoAYBF1FVDbOtv1xjnXPz5HVwkASLq6Cqitb12r1c0r6CoBAClQVwG1orFBP3ZTux4/fFznz7OIIQAkWV0FlBR0lXjh5TP62ndOV7sUAMAl1F1A/ejb1qnBpINPf7/apQAALqGqAWVmwwW3+82sz8wGl/M9W1Zeoe63ttJVAgASLlJAmdkWM7vLzD5kZveGt7cs8TX6JGXC2/2S5O4TknLhY8tm2+Z2PfXdF/W9068t59sAACqwotwdzWyTpCFJmyRlw5+cJJPUIelOM8tIOiJp2N2fv8RrZcLn5/VIeii8nZXUJWliwWJeeE760/devO0dPyfduks686r06R3zn7PlF6Rb3i+9clIfOPyr2nLFadmf/YG0qjms4FekH75bOn1MenT3/Off9uvSTT8ZvPff/Ob8x//Nh6SOO6R/+Zr0uQfmP977Een6d0n//GXpwH+Z//hP7JHe/E7pyOPSFx6c//hP/6F07Y3S4b+TvvhH8x+/a6+0ZoP09Uekyf85//F7/ly6uk069Gnpyf89//H375euWCl9ZZ/01GfmP/7Lfxv8/qePS88+dvFjTc3SLz4S3P78iJT9/MWPr1wr7fyL4PbE70vfnrz48dXXSXfvC27/3f3S9/7fxY+3dUg/8/Hg9l//hnTyyMWPv+lfST/50eD2I7ukF7978eNv6ZH6fj+4/dAvSq8WdbXP/Kj0o+GB+1/cLb1R9D8ub/tx6fbfCG4X/7uTlvRvTw9/YP7j/Nvj355Un//2FlFWQJnZLklrJA25+yVnF5jZGkkDZnbK3T+1wG4Zd58ws/z9lqLH20q87oCkAUl653VXlVP2gq5qatSVKxqUe/WM1ucDCgCQKOZ+6enWYThNuPvRJb1wcMTVWxxSZtYXnsqTme139x3hWNR4GFp9kra7+9BCr93d3e1TU1NLKWee3//rp/SXk/+sJz9yp5qbGit6LQBA+cxs2t27F9tv0TEod9+31HAKn3d0gSOo2XAyRL+kjJl1SZrUhaOojKTxpb7fUm3rbNdrb5zXl46cXO63AgBEEHWSxJ6ob+juM+ERVKvCUHL3MQVh1RfeX3j8KSbvyrRq5RWNOvAM080BIInKniRRZLeZHZH0sLu/GOUF3H1U0mjB/ZGItURy5YpGvefGa3Xw6ePyn3UVjIcBABIg6nVQu8LTdz3hFPON8ZV0+fR2rtd3T7+mZ773UrVLAQAUiRRQ7v5I+PuAuz+q4PTcY2Z2V6zVLbMf61wniUUMASCJIl+om/9tZg9LGpN0VNJRM7vbzLbFWOOyaV/VrJs3rNEB2h4BQOJEHYMaM7P8/PQRd7+n4LFDkmRm917iOqjE2Na5Xn944FmdfPl1tV1zZbXLAQCEoo5B5STd5+43uvu+4gfDU32zFVV2mfRubpe79MRhFjEEgCSJGlB73P3AJR6/Uxe3Mkqsd1y3WutXX8k4FAAkTEWTJPKKZ/G5+33u/mT0si4fM9O2znZ94dkTOnP2fLXLAQCEok6SuHf+JustsT0VtnWu10uvn9XU86k4KwkAdSHqKb6LmruGbY0udcov0W6/oU1XrGhgjSgASJClLLexS9JWSWsldZlZT9Eu+SU0Ej9zr9jKK1boto42HXzmuH73p95e7XIAAFpCQIWz9faFK96apIeLdsm6+6E4i7ucejvb9buffUrZEy8rs+6aapcDAHVvyddBufuImfWm+ZReKXd0tkuffUoHnzlOQAFAAkSdxVcynNI6SUKSNqxdqc43rdKBpxmHAoAkKHdF3T+RtN/dD4b3Hyu1m4IxqtSNQeVt62zX6BeyOv2DN7TmqqZqlwMAda3cU3zFa1G0SSpe8dYkDVZcURX1bm7XJ544on947oR+6p3XVbscAKhrZQWUu99XtGlXqQkRZpbq5Wm3vGWt1q5s0sGnjxNQAFBlUcegFpqt11tBLVXX2GC646Z2PX74uM6d98WfAABYNuWOQZUac5q3m4IxqAcrqqjKtm1u16OHvqMnv31KW9/aWu1yAKBulTsGVWrMqVjqx6Ak6T03rtOKBtOBp48TUABQReUGVMkxp2JpH4OSpDVXNalnY6sOPnNcgz/RWe1yAKBulTUGtVg4mdkaM7tbUk0M3PRubtcz33tJx069Wu1SAKBuRW0WK0kys9XhUhtrJU1L2hlDTVW3rbNdkvQ4zWMBoGqiLrfRa2azko4qCKaZ8PdkjLVVTWbdNdp07dV0NweAKop6BNXn7q3u3iZpwN1bFXQzz8VXWnVt62zXF4+c1Ktnzla7FACoS1EDaqrg9lpJcvfTlZeTHHe+fb3OnD2v//jQVwkpAKiCyGNQZnZXeHOtmd0c3u6qvKRkeFemTb/z3s167Bvf0z17v6TvnX6t2iUBQF2JGlBZSb8dTpAYlfRIOMW8I6a6EuHe92S075e6dfTEK/rZP/5Hff07NXWQCACJFrnVkbt3u/vz7n7a3W9QMC71wZjrq7q+t6/X2Adv04qGBu345Jf0ua//S7VLAoC6UNE080Lufig8oqo5m9+8Wn/1a7fppjet0n1/MaNPPPFNudfEJV8AkFgVB1R4LdRqM1utxdshpVb7qmb95cC79dM3X6eRzx3Wb+3/ql4/e67aZQFAzYp6HdQuMztvZucUTC0/Ff4eiLO4pGluatTH37dFv9l3ox6d+Y5+8VNf1uwrZ6pdFgDUpKhHUB2S1rp7o7s35H9L+liMtSWSmek3+96mj//bW/S1Y6f1s3/8j3ru+y9VuywAqDlRA2p8geue9lRSTJr8zM3X6S8H3q0fnDmvuz7xRX3+2RPVLgkAakrUgPIFJkTsil5K+txy/Vp99tdv1w+tvUq/8r8m9edfer7aJQFAzYgaUPdJGjez58zssfBnStIDMdaWCj/UcpXGPnib7rhpnT7y2af0e5/9us6eO1/tsgAg9cpdD6pYRtL9urj3Xk0sWBjFNVeu0N5f6tZH/+5p7fuHozp68lX90S/cotXNTdUuDQBSK2pADbn7geKNtbBgYVSNDab/9N63q2PdNfqdz3xdd3/ii/of/65H17etrHZpAJBKUTtJzAun0NYKaqkJ77v1ev35v79Vx196XT/3iX/S5POz1S4JAFKprIAysz8xs20F9x8r8fP3koaXrdIUua3jWn3m125Xy1VNev++L+uR6WPVLgkAUqfcU3xWdL9N87tG1O0YVCmbrr1af/Wrt+uDn57Wb+3/qo6ceFkfuvMmNTQU/6cEAJRSVkC5+31Fm3a5+6Hi/ep5DKqUNSub9Ge/cqs+8tmn9Iknjih74hX9wc6btfKKqEN/AFA/InczX+ChIxXUUpOaGhv0X3/+h1lbCgCWKLZu5qF9Mb9eTTAz3fuejD71AdaWAoBylTtJ4ryZnVvsR1L/Mtebar2bWVsKAMpV7hHUWNgQttHdGyXdKemGom094XZcwuY3r9Znfu12db6ZtaUA4FLKHa0v7rG3xt2PFm5w95nCqehY2LpVV+r/7Hq3Bse+ppHPHdY3j7+s+3+ikxl+AFCg3Fl8xQMmrQvs2lJZOfWjualR//19W9Sx7hr9t4ln9ejMd6pdEgAkStT5zjeY2RZ3fzK/wcy2KDjN92gsldUBM9N/6LtR78q0sqYUgLrxgTJbOkQKKHe/38z+3sw2KWgYm5GUldQb5fXq3bszbXp3pq3aZQDAZfGBMveLfMWou99pZrdI6paUvUR/PgAAlixSQJnZXQpC6ZCkhS7aBQAgsqgX6r6v1EYzW11BLQAAzIkaUA8pGHMqNlBBLQAAzIk6BrVd0kfNLKsLq+qagkkSD8ZRGACgvkUNqG5JI5KKV+PjOigAQCziXvI9V2pnAACWKup1UBeFk5mtkdQnltsAAMSkouU2zGy1mW2UtFbStKSdMdQEAEDk66B6Je2XlG/DbeHt4qayAABEEvUIqs/dW929TdKAu7cqaHfEGBQAIBZRA2qq4PZaqWTHcwAAIos8BhW2O5KktWZ2c3i7q/KSAACIHlBZSb8dTpAYlfSImZ2U1BFTXQCAOhd1mvkhBRfr5t1gZreE2wEAqFglp/i2FG3ykjsCABBBpIAys7sljRWF1Gkz2xZPWQCAehe11VHO3W8o3ODuR8MFDBdlZn3hze3uPhRu61cwTb3L3Uci1gUAqBFRT/GtWWB762JPNLMuBcE0IanLzDJhOCnclisIMABAnYoaULcWj0GFp/e2LvZEd59x9yEza1GwKm9WUo8urC+VFdPVAaDuRT3Ft0fSATPbpCBQMuHv3iW8RrcudJ4oXqajrXhnMxtQuCDi9ddfv9R6AQApE+kIyt1Pu3u3guawo5Lucfced39xCa8xIamlYOzpkqcH3X3U3bvdvXvdunVRygYApEjUIyhJcyGzJGY2LOmIu4/qQjBN6sJRVEbSeCV1AQDSb9EjKDPbFXaMWBIz22Rm95Z4aK+kbDgRoiU8MhqTlMlPjogSfACA2rLoEZS77wtDKiNpr7s/f6n9zWy1pN+W9BV3/1SJ18vqwoSIiYLtTC0HAMwp6xRfGFKbJN0XXuuUVXB67oiCU3Nt4e+OcNuIux9dnpIBAPWg7DGoMHDul4LTdwrGijKSTks6qmDKOL34AACxiNos9qiCUDoQbzkAAAQiN4sFAGA5EVAAgEQioAAAiURAAQASKdIkifDC3RZ3fzK87mlAktz9wfhKAwDUs6hHUPcrmGIuBTP52hQ0j/1QLFUBAOpe1F584+7+aHg91FZ375GksNsEAAAVi3oEdSr83SdprGC7V1YOAACBqEdQW81sraQhheNPZtarMlbUBQCgHFHXg/qYgjDa7e4Hw3BiFVwAQGwirwcVNpBdHS79PunutD0CAMQm8nVQZvZJBR3ND0o6ZWYPxVYVAKDuRQooM/uwgpl8De7e6u6Nkh5mmjkAIC5Rj6Cy7v5I4Ybw/unKSwIAIHpALTSd/GTUQgAAKBQ1oDrCFkdzwvZHt1ZaEAAAUvRZfKOSDpqZS5pVMOW8RdLWuAoDANS3qCvqnpbUbWZ3K+jJN29MCgCASkTtZn6XCCUAwDKKOgb1vlIbi8elAACIKmpAPSQpW2L7QAW1AAAwJ+okie2SPmpmWQXdJCTJJPVKYtFCAEDFogZUt6QRBTP4CrVUVg4AAIGoATVUqjmsmeVK7QwAwFJFXW7jQNjJ/N7wAl2Z2RY6mgMA4hK1WWyvgi7mWxVcByVJp81sW1yFAQDqW+RJEu7eLc2Fldz9qJmxaCEAIBZRp5l/ZYHtCzWRBQBgSaIG1K1mtiq87RLNYgEA8Yp6im+PpENmdkqSzKxFwfVQvXEVBgCob5U0i72BZrEAgOUS9QhK0twqujKzW8xstbu/GE9ZAIB6F3Wa+Ukz+/mCTTlJO82MFXUBALGIegT1KUkfM7Pt7v6r7n5U0r5wLAoAgIpFncX3ZXe/QVKjmT1nZjeH20t1OAcAYMmiBlSrJLn7bkkflPS4mX1IXAcFAIhJ1IC6L9+Dz90nFMzk+3FJ++IpCwBQ76IG1C5JHWa2RZLcPefu2yXdH1tlAIC6FvU6qEP522a2zd0Phts5ggIAxCLqEVShvTG8BgAAF4kjoCyG1wAA4CJxBBQz9wAAsYsjoI7G8BoAAFyk4oBy9zvzt/NTzwEAqFTFAWVmq/M/koZiqAkAgMjNYneZ2XkzO6egUeyp8PdAnMUBAOpX1COoDklr3b3R3RvyvyV9LMbaAAB1LGpAjYeLFhbbU0kxAADkRQ0oX2BCxK7opQAAcEHU9aDuk3SLmUkXlthok7RJ0oMx1AUAqHNRAyqjoDFsrmCbSRqsuCIAABQ9oIbc/UDxRpZ8BwDEJdIYVD6cwuuftoTXQF3U5RwAgEpEvlDXzD6p4BTfQUmnzOyh2KoCANS9qBfqfljBVPMGd29190ZJD4fLvgMAULGoR1BZd3+kcEN4v9S1UQAALFnk66AW2M4kCQBALCK3OspPjMgLL9y9tdKCAACQok8zH5V00Mxc0qykVkktkrbGVRgAoL5FCqiwD1+3md2t4KLdeWNSAABUIuoRlKS5iRFzzGyjuz9fUUUAAIgFCwEACcWChQCARIp6ii+/YOFF1z2Z2UcXe6KZtUjqC+/2uPtQuL1fQch1uftIxLoAADWiGgsW3iOp1d3HJMnMBsJwkrtPSMqZWd+lXgAAUPsu+4KF7j7q7qPh3YykCUk9urCuVFZSV8S6AAA1omoLFppZRtKsu2fD036F2krsP6BwjOv666+PWDYAIC2quWBhv7vvDm/nFFzsu6DwqGtUkrq7uxdqtQQAqBFVWbDQzPrzEyHC8aZJBZ0opCD8xiPWBQCoERUtWFhi+6ILFoaBNGxm02Y2HT5vTFImPzkinCwBAKhjFXWSiCIMn44S25laDgCYU3EnCQAAlgMBBQBIpLICysw2LnDdEwAAy6LcI6hhFVw8a2ZbSu1UvIghAABRlRtQR9390YL7C7UieqDCegAAkFT+LL5JM/umpCPh/YyZbS/axxSsqEtIAQAqVlZAufsjZjYhqTvctEPS/qLdltpJAgCABZV9HVTYvfyAJJnZbKmLcsvtJAEAwGKidpI4JM2tprslPzminE4SAACUI/J1UGb2SQVNXg9KOmVmD8VWFQCg7kVd8v3DChYtbHD3VndvlPSwmX0o3vIAAPUq6hFU1t0fKdwQ3i+1yi4AAEsWeUXdBbYzSQIAEIuoAdVR3DUibIV0a6UFAQAgRV9uY1TSQTNzSbMKVsNtUXChLgAAFYsUUOE1Ud1mdreCFXDnjUkBAFCJihYsJJQAAMuF9aAAAIlEQAEAEomAAgAkEgEFAEikqK2O7lpoVV0AAOIQ9QjqfaU2suQ7ACAuUQPqIUnZEtsHKqgFAIA5Ua+D2i7po2aWVbDkhhSsqNsr6cE4CgMA1LeoAdUtaURBm6NCLZWVAwBAIGpADbn7geKNZpYrtTMAAEsVdcn3A+Fy7/eGXcxlZltKhRYAAFFEnWbeq2Cp960KmsVK0mkz2xZXYQCA+hZ5koS7d0tzYSV3P2pmXbFVBgCoa1GnmX9lge0LrbQLAMCSRA2oW81sVXjbJVbUBQDEK+opvj2SDpnZKUkysxYF10P1xlUYAKC+VbKi7g2sqAsAWC6sqAsASCSW2wAAJFLkgDKzXWY2a2bnzOw5M7s3zsIAAPUt0ik+M9ul4CLdHQq6mq+VNGBms+7+aIz1AQDqVOQxKHe/r+DuUUn3mdmHKy8JAIDop/iKu5jnlVojCgCAJYsaUHSMAAAsq7JO8ZnZQ0Wb1prZsKSZgm0ZSbviKgwAUN/KHYPqkLRXC5/ay2utrBwAAALlBlTJBQqLsdwGACAuZQXUAqvn3qX5R0y7JfXEUBcAoM5FvQ7qkwrGnHK6cNqvVVJLTHUBAOpc1Ougpouug5I0dwEvAAAVi/s6qPGohQAAUChqQM2Y2TYz22hmq/M/kobiLA4AUL+inuLrUzDtPM8lWfj7g5UWBQBA1COoFklr3b0h/Gl09wZJH4uxNgBAHYt8ii9cVbfYnkqKAQAgL3IvPjPbWGI7s/gAALGIOgZ1n6RbzEy60MG8TdImSQ/GUBcAoM5FDaiMpPsVXKibZ5IGK64IAABFD6iSvfnM7GSF9QAAICniGNQlGseeqqAWAADmRAqowotziy7UHY65PgBAnYp6ii+nCxfn5rkuHpMCACCyqNPMx/IX5+Z/JHUr6DABAEDFogbUvOud3P2QgmnmAABULOokiVJdJCSWfAcAxCTqgoWPldickTRWWTkAAASiTpJoU9B3r3BSRNbdj1ZeEgAA8V+ou9Hdn6+sJAAAYrhQlwULAQDLIeqFurvM7LyZnVNwmu9U+HugzOf3mdl40bb+cDv9/AAAkaeZdyhYsLAxyoKF7j5ReN/M+gu258yM66kAoM5FDajxmBcs7NGFZTuykroivg4AoEYkZcHClqL7bcU7mNmAmU2Z2dSJEycivg0AIC2SsmBhTotc5Ovuo5JGJam7u9sjvAcAIEWSsmDhpC4cRWUkjV9iXwBAHajKgoXhpIhuM+t39zF3HzOzwfzkiOJJFACA+hMpoBZasDBsGFvO88dU1BbJ3Uei1AIAqE1RJ0kAALCsCCgAQCIRUACARCKgAACJREABABKJgAIAJBIBBQBIJAIKAJBIBBQAIJEIKABAIhFQAIBEIqAAAIlEQAEAEomAAgAkEgEFAEgkAgoAkEgEFAAgkQgoAEAiEVAAgEQioAAAiURAAQASiYACACQSAQUASCQCCgCQSAQUACCRCCgAQCIRUACARCKgAACJREABABKJgAIAJBIBBQBIJAIKAJBIBBQAIJFWVLsAACjHG2+8oWPHjum1116rdikoU3NzszZs2KCmpqZIzyegAKTCsWPHtGrVKm3cuFFmVu1ysAh318mTJ3Xs2DFt2rQp0mtwig9AKrz22mtqa2sjnFLCzNTW1lbRES8BBSA1CKd0qfTvRUABABKJgAKAMkxMTGj79u1Ve36xbDarHTt2LPj42NhY2fsmFZMkAKTOf/6bp/SN774Y62u+/brV+r2ffseCj/f19Wnv3r2RX7/S5xfLZDLav39/ycdyuZzGx8fV39+/6L5JxhEUAEQwOjqqmZkZjY6Ozm0bGRnR2NjY3O9SJiYmNDIyopmZGUlBmIyMjGhiYkKjo6MXHWmNjIxoaGho7nnbt2+fe/2ZmZm5/WZmZjQxMTH3k81mNTU1pYmJibnHC4/e8s9fqMbEcPfU/WzdutUB1JdvfOMb1S7B+/v73d19eHjYp6en3d19fHzc9+7d6/v37/f9+/dftF+xrq4ud3c/cuSIDw4Ourv74ODg3GsNDAxc9PzC/dzdM5lMyXoGBwd9fHx87jmlasjfL6fOOJX6u0ma8jK+6zmCAoAlGh8fV0tLi6Tg9Nn4+Li6uro0OTmpbDa74FhTJpOZt21mZkazs7OamZnR7t27L/m+XV1dJbc/8MADGh8fV0dHh3K53KK15+tI+mk/AgoAlqirq0vZbFZSMAGhp6dHkrRz505J0sDAQNmvlQ+zrq6ueQGWf4/FTExMaHh4WNPT03On9fLypxLzOjo6NDs7K0mLhlm1MUkCAMowMzOjmZkZZbNZDQ8Pa2T6w8rHAAAHdUlEQVRkZG774OCgstms9uzZo9bWVs3OzuqBBx646Iin8PkTExOamZlRLpfT4ODg3GtJwWSKnp4eTUxMKJfLzf2empqae42urq6LXm9ycnLu+YUTI8bGxtTX13fRvoODg3PjWrlcbm7/JLLgdGC6dHd3+9TUVLXLAHAZPf3009q8eXO1y1jQ0NCQhoeHF7xfr0r93cxs2t27F3suR1AAEIOdO3dqbGxMmUxG2Wx27nQfoiOgACAGXV1dc6f0FprMgKVhkgQAIJEIKABAIhFQAIBEIqAAAIlEQAFAGYq7kcfdnTyKWu+Qziw+AOn0p++dv+0dPyfduks686r06RJfnlt+Qbrl/dIrJ6WHP3DxY7/8t5d8u+Ju5H19fVW/zqnWO6RzBAUAMcsf2RR3Lpfmd0Ev7lJe/NyxsTFNTExo9+7dC7YmqtkO6eV0lE3aD93MgfqTpG7meX19fQvuW6pzeaku6O7zu5QXPjf/nnv37p3rWL7Y+ySpQzrdzAGgyoaGhjQ0NDR3ZFSqc3mpLujS/At7C5+bv93a2lryfWu5QzoBBQAxGB4e1vDw8CU7mS/UBT1utdIhnYACgDIUdgSXghltU1NTJcdcFupcPjw8PDe2k++Cnr+d/9Iv9dxsNqvx8fG5I67F3mdwcHDuffKNtUt1SL/Ue09OTs6NfxV3SM9kMvM6pI+Pj8+NV8WFbuYAUiHp3cxRWiXdzDmCAgAkEgEFAEgkAgpAaqRxSKKeVfr3IqAApEJzc7NOnjxJSKWEu+vkyZNqbm6O/Bq0OgKQChs2bNCxY8d04sSJapeCMjU3N2vDhg2Rn09AAUiFpqYmbdq0qdpl4DJKTECZWb+knKQudx+pdj0AgOpKxBhUGE5y9wlJOTPrq3JJAIAqS0RASeqRlO+/kZVUujkUAKBuJOUUX0vR/bbiHcxsQFK+ydXLZnZ42asKXCvphcv0XpdDrX0eic+UBrX2eSQ+UyXeWs5OSQmonKTSrXpD7j4qafTylHOBmU2V05IjLWrt80h8pjSotc8j8Zkuh6Sc4pvUhaOojKTxS+wLAKgDiQgodx+TlMlPjggnSwAA6lhSTvEpwVPLL/tpxWVWa59H4jOlQa19HonPtOxSudwGAKD2JeIUHwAAxQgoAEAiEVAlmFmLmfWHP8PVridutfKZzKwr/3eqdi1xCT9PX3jdX2qFn2G8aFv+sw1Wq65KFH+mWvieKPV3Knis6p+JgCrtHkmt4exCpf3LolA4UzJT7TpisrtgBmjqP1P4t8mGs1izZpbajirFM3FroZ1ZidnFqf+eWGjGdFK+JxIziy9JwouC8zKS9larljiFX+LZRXdMgfDLYNrMMgmeAbpUUwo+0w5JmRq73KJH0kPh7Xw7s1R/Pr4nlh9HUJcQ/qFm3T0Rf6wYZGros3SEP7NmttfMittlpY675xR8ye1X8NlqyaLtzNKK74nlQ0BdWr+77652EXEws74a+z9ySToSfqlP60KfxtQKT4NNuHtHwf1asWg7sxTje2KZEFALMLP+/KmjNJ4vL2E2HBDtVzBmk9rxjdBkwe0WBV+AaZdx95nw9h7V1hd6TbYz43tieRFQJYT/0IbNbNrMpqtdTxzcfSb8P6NWzT/dkjrhwHRLQXusRF0BH9GomQ2En+meNH+m8Auuu2ByROrbmRV/plr4nijxd0rU9wSdJAAAicQRFAAgkQgoAEAiEVAAgEQioAAAiURAAQASiYACACQSAQUASCQCCihD2O8vUc1Ak1gTECcCCijP/vBnzuVcXmGB95pXE1BLWG4DKMMCrXm2XsYS5r1XGtsFAUvBERQQQbgq7GVp5no53wtIEo6ggEWEHZ33KVjzZ3vYJLRHUlcYHrl8Y9dwXaoHFHTv7pE07u4T+caiChYlHJe0U9JD7j4Wvn5+9dK554SvV/K9imsqeO8BXVhsLlPUaXs4fCw/brVdwZIlqW1Ki9pGs1igDOEX/FBBGPRL2unuO4r2OyJpa7hOlcIu173ungvHkYYUnK7LSEH36HCfPQVLh59y97UFr7nQexXXNPdeBc/bnl+rKLw/HG7LhoF2yt0txv9UQGw4xQfEJAyAXD4gQlOS7glvz0rKunsuXNYgv/bTjnw4FbzWkpY6KFguYe69w9ccKHitXLg9W7wvkESc4gPik5HmLVy3XxdOuUmlF1aczZ++UxBiUjDmtJQAyRQ8t1AufCwfhvOW8jazFsIKSURAARUys/z4UVZSS4TZddMKjqJmwtfbF/6eFxz598ofBRXIqvREihaVCCUgDTjFB0ST1YWJDRl3z4an1HIFgSUzu+Sy2eFjrQXhVHhqL38kNu+9il+nYIXhwvfulzRWFHLMBkRqEFDAIsIQGVKwNPaAFExukJQN72cKdu+VtNvM+sOAyIQTIfok7ZbUZ2aD+SAKX+fhcFufpG5Ju8J9swu9V6maFEy+yL/3QPjeO4r2z+Tf38yGw+cNFwYbkBTM4gMAJBJHUACARCKgAACJREABABKJgAIAJBIBBQBIJAIKAJBIBBQAIJEIKABAIhFQAIBE+v8GuPTiAB4I/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Show plot inline in ipython.\n",
    "%matplotlib inline\n",
    "\n",
    "# Plot properties.\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')\n",
    "plt.figure(figsize=(6,6))\n",
    "\n",
    "# Plot the two data series.\n",
    "plt.plot(range(1,1+NUM_RUNS), nnzs_log, label='log heuristic')\n",
    "plt.plot((1, NUM_RUNS), (nnz_l1, nnz_l1), linestyle='--', label='l1-norm heuristic')\n",
    "\n",
    "# Format and show plot.\n",
    "plt.xlabel('iteration', fontsize=16)\n",
    "plt.ylabel('number of non-zeros (cardinality)', fontsize=16)\n",
    "plt.ylim(0,n)\n",
    "plt.xlim(1,NUM_RUNS)\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
